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Abstract 

For the first time exact analytical solutions to the eikonal equations in (1+1) 
"^ ■ dimensions with a refractive index being a saturated function of intensity are 

constructed. It is demonstrated that the solutions exhibit collapse; an explicit 
analytical expression for the self-focusing position, where the intensity tends to 
2 ' infinity, is found. Based on an approximated Lie symmetry group, solutions to 

the eikonal equations with arbitrary nonlinear refractive index are constructed. 
Comparison between exact and approximate solutions is presented. Approximate 
solutions to the nonlinear Schrodinger equation in (1+2) dimensions with arbi- 
trary refractive index and initial intensity distribution are obtained. A particular 
^SJ ! case of refractive index consisting of Kerr refraction and multiphoton ionization 

C^ ■ is considered. It is demonstrated that the beam collapse can take place not only 

at the beam axis but also in an off-axis ring region around it. An analytical 
condition distinguishing these two cases is obtained and explicit formula for the 
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1 Introduction 

The Lie symmetry analysis of differential equations finds a great number of ap- 
plications in mathematical modeling of physical problems nowadays (see, e.g. 
Refs. [IH2])- )• Nonlinear optics certainly occupies a particular place among 
these. 

In 1960-ies it became evident that for an adequate mathematical description 
of the process of highly intense light propagation the refractive index has to de- 
pend on the intensity of applied electric field n = n{I). For moderate intensities 
achievable at that time it was sufficient to use the so-called Kerr form of the 
refractive index n = hq + 71-2/ with uq being of the order of unity and 71-2 vary- 
ing from 10^^^ to 10^^^ in (W/cm^)^^ depending on the material. The basic 
mathematical model for the propagation of intense monochromatic light that is 
successfully applied for a long time (see, e.g. classical monographs [5| Chapt.17] 



and [5], or [SHS] for recent achievements) is the nonlinear Schrodinger equation 
(NLSE) or its approximation in the hmit of geometrical optics, the eikonal equa- 
tion. 

For the first time, exact analytical solutions to the eikonal equations with 
Kerr-type refractive index in (1+1) and (1+2) dimensions were constructed by 
Akhmanov et al. in Ref. [5l[10]. The authors demonstrated that in both cases 
the solutions exhibit singularities at certain points and found explicit analytical 
expressions for them. Later in Ref. |lllll2j it was demonstrated that these solu- 
tions can be derived in a regular manner using the Lie-Backlund symmetry group 
admitted by the eikonal equation with Kerr refractive index. 

Lie symmetry group analysis of NLSE has been performed by many authors 
(see, e.g. [H Chap. 16]). In particular, L. Gagnon and P. Winternitz in Ref. 
[T3HTH] found exact solutions in (1+2) dimensions, however, these solution did 
not correspond to localized (symmetric) intensity distributions typical for usual 
experimental conditions. A set of two coupled NLSE was analyzed by means of 
Lie group technique and the general Lie group of point symmetries, its Lie algebra, 
and a group of adjoint representations that corresponds to the Lie algebra were 
identified in Ref [16]. 

For the initial conditions actual for typical experiments, analytical solutions to 
the eikonal equation in (1+1) dimensions based on the symmetry group approach 
were obtained in Ref. |12J. In this paper approximate solutions for various initial 
intensity distributions and Kerr-type media were constructed. Later, using the 
similar group analysis technique approximate solutions in (1+2) dimensions for 
arbitrary initial intensity profile and the same form of the refractive index were 
found in Refs. [I7l[l8]. Based of the obtained solutions, the authors proceeded 
to investigate a global behavior of the solutions and to get explicit analytical 
expressions for the nonlinear self-focusing position and the value of critical power 
required for beam collapse. 

On the other hand, modern experimental facilities allow one to achieve very 
intense laser beams leading to highly nonlinear media response. In such situa- 
tions the Kerr approximation to the refractive index ceases to be sufficient and 
higher order terms with respect to power of the light intensity must be taken into 
account. 

In the present paper we study the problem of light propagation in media with 
highly nonlinear response. Based on Lie symmetry group analysis, we constructed 
an exact solution to the eikonal equation in (1+1) dimensions for a special higher- 
order form of the refractive index, and also approximate analytical solutions to 
the problem in both (1+1) and (1+2) dimensions with the refractive index being 
an arbitrary function of the intensity. 

The paper is organized as follows. First, after describing the model equa- 
tions, we consider the problem of light propagation in (1+1) dimensions under 
an approximation of geometrical optics. Based of the formalism of Ref. [llj we 
construct the Lie-Backlund symmetry group admitted by the eikonal equation 
and consider such a superposition of the symmetry operators that yields a local- 



ized initial light intensity distribution. The use of this combination of operators 
gives us an exact solution of the eikonal equation with the refractive index being 
a saturated function of the intensity of the applied electric field. The solution 
exhibits a singularity: the on-axial intensity asymptotically tends to infinity at a 
certain propagation distance. 

In the next section we construct an approximate analytical solution to the 
eikonal equation in (1+1) dimensions with arbitrary nonlinear media response. 
This solution is obtained on the basis of the approximate Lie-Backlund symmetry 
group. In order to test the applicability of used approximation the solutions 
obtained on the basis of exact and approximate Lie symmetry groups under the 
same initial conditions are compared. 

The last section is devoted to the construction of approximate analytical so- 
lutions of the Schrodinger equation in (1+2) dimensions when both the refractive 
index and the initial intensity distributions are arbitrary. The final result is a 
system of algebraic equations which has to be resolved for every particular initial 
condition and nonlinear media response. The most typical situation in modern 
experiments (see e.g. [1119]) is the propagation of the Gaussian beam in an ion- 
izing media, with the refractive index being a polynomial function of the light 
intensity. Therefore, we consider this problem as a particular application of the 
obtained results. Influence of the higher order nonlinear term in the refractive 
index on the beam collapse is considered and result is compared with the previous 
one obtained in Refs. [TTlfTS]. 

2 Model equations 

Let us start from the NLSE: 

i£, + ^V\£ + kQn{\£\'')£ = ^. (2.1) 

2^0 

Here £ is the slowly-varying envelope of the electric field, z is the propagation 
length, ko is the wave number k^ = uqujq/c, wq is the carrier frequency of the laser 
irradiation, c is the velocity of light and n = n(\£\^) is a nonlinear refractive index 
in a general form (see e.g. [HIEIIH])- Due to high intensities of light available in 
modern experiments the refractive index becomes highly nonlinear n = n(|i?p). 
The Laplace operator V^ is usually responsible for light diffraction. Explicitly, 
^\ = dxx or dxx + dyy for (1+1) or (1+2) dimensional cases correspondingly. 

Let us now represent electric field £ in the eikonal form: £ = \/l eyip^ikoQ). 
Then, starting from Eq. (12. ip . after some algebraic manipulations we obtain 



Qz = -liQ.? + n{I) + -^2 f^dxix'^-^dx^) ] , (2.2) 



-dx{IQx)-(.i^-l)^-^, (2.3) 



where v = \ and u = 2 correspond to the (1+1) and (1+2) dimensional cases, 
and X denotes the transverse spatial variable. 

Let us differentiate the first equation with respect to x and introduce a new 
variable v = Qx- For the sake of convenience, we introduce diniensionless vari- 
ables / = ///q, X = x/ii^in, z = z/tt^in, where /q is an initial peak intensity of the 
light beam and Wm is an initial beam radius. In what follows the dimensionless 
parameters shall always be used, omitting the tilde for simplicity. Moreover, let 
us introduce new dimensionless variables a = n2lo, = (2^o^fn)~^ in diffractive 
case. Thus, finally we get the following equations: 



Vz + vvx - aiplx - Odx 



A-u 



^/7 
vl 



a. (x^-'d^Vi' 



0, 



Iz + via: + IVa: + {u - 1) =0, (p = dlU . 



(2.4) 
(2.5) 



Evidently, Eqs. (|2.4p . (|2.5|) must be supplemented with a boundary conditions. 
In case of collimated beam these read 



v{0,x) =0, 



/(0,x)=/o(x), 



(2.6) 



In several cases the term with higher order derivatives can also be neglected 
Refs. [T9j and equations (|2.4p . (j2.5p acquire a rather simple form: 



Vz + vvx — acplx = , 



Iv 

Iz + Vix + IVx + (Z^ - 1) 

X 



0. 



(2.7) 



Further simplification of these equations can be performed in (1+1) dimensions 
if one notices that in this case the system (12. 7p is linear with respect to the first 
order derivatives. Therefore, it is convenient to use the hodograph transformation 
in order to transform it into a linear system of partial differential equations. In 
doing so, in (1+1) dimensions one obtains 



OT,, 



V{I) 



X/ = 0, Xv + ri = Q. 



In (1+z^) dimensions the Eqs. ()2.7p read 



aT„ 



V{I) 



XV+T-I + 



{v-l)v 

Xl + TV 



Xv (Iti - t) + TTI 



T 

T 
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I2 



^{I) 



The boundary conditions are transformed as follows: for v = 

r = 0, X = H{I), 



0. 



(21 



(2.9) 



(2.10) 



where H(I) is a function inverse to a smooth initial intensity distribution Iq 



Io{x). Evidently, for example, in case of a Gaussian beam we have x = Y^ln(l//). 



3 An exact solution to the eikonal equations in (1+1) 
dimensions 

Let us now construct a new exact analytical solution to Eqs. (j2.8p . In the present 
paper, for the first time, we consider a special form of the nonlinear media re- 
sponse with the refractive index being a saturated function of intensity. Let us 
rewrite Eqs. (|2.8|) as follows: 



i^{i)xi = 0, 



Xv + Tl = 0, 



(3.11) 



where ijj = I/{a(p). As a boundary condition, we take a collimated continuous 
wave beam with a localized symmetric intensity distribution at the entry plane 
of a nonlinear media 



t(/,0) = 0. 



xiI,0) = XoiI). 



(3.12) 



Our goal is now to construct an exact analytical solution to the system of 
equations (jS.lip with a nonlinear function ip{I) that corresponds to a saturating 
dependence of the refractive index on the intensity. For this goal, ip is taken in 
the form e ' ° /n2lo. The refractive index corresponding to this choice of tp is 
presented in Fig. [TJ 




Figure 1: Nonlinear refractive index corresponding to "0 = e /a as a function of 
the normalized intensity /. 

Let us first sketch the broad outlines of our solution: first, we construct a 
Lie symmetry group admitted by Eqs. (j3.1ip . Second, following a formal scheme 
reported in Refs. [20], the obtained group shall be restricted to the surface of 
boundary conditions: v = 0, t = 0. All derivatives of r with respect to / have to 
vanish too. Third, based on the requirement of vanishing of canonical coordinates 



of the group generators on the boundary, we shall construct such a linear super- 
position of them, which provides a localized beam intensity distribution. Finally, 
the integration of the constructed superposition shall yield a desired solution to 
Eqs. (I3TT]1 . 

3.1 Recursion operators and Lie-Backlund symmetries of the 
second order 

We start from the search for the Lie-Backlund symmetry group admissible by 
Eqs. (j3.1ip . It is generated by the canonical infinitesimal operators [T] 

X = fdr + 9''d^, (3.13) 

with coordinates f and g'^ . For the Lie-Backlund symmetry of arbitrary order 
s > 1 the coordinates /^ and g^ depend on v, /, r, x and corresponding derivatives 
of r and x up to the s-th order with respect to / 

/' = f'{v, I, T, X, -, r/, xf ), g' = g'iv, I, T, X, -, r!,x'i) ■ 

Here, the index s stands for the order of the derivatives: rf = d'^r/dP , etc. The 
coordinates f^ and g'^ are found from the determining equations that in the case 
of Eqs. (I31T]) read [II]: 

D,if)-i;Di{g')=0, D,{g') + Di{f) = , (3.14) 

where Dj and D^j are operators of the total differentiation with respect to / and 
v. 

oo 

Z' (3-15) 

In order to solve the Eqs. (|3.14|) . it appears more convenient to use a recursion 
operator [H]. The latter is defined as 2 x 2 matrix operator transforming any 
linear solution of the determining equation (I3.14p of the order s to the solution 
of these equations of higher order (s + 1) 

Substitution of Eq. ()3.16p into the determining equation (I3.14p yields the 
following system of equations for the elements of L 



{DjL'' + D,L'')f + {DiL'^ + D,L^^)g' = 0, 
which should be valid for any solutions /^ and ^f*. 



(3.17) 



Explicit formulae for the components U^ {i,j = 1...3) of recursion operators 
()3.16p are given in Ref. [11] 



Li = ( ;: ^J'" 1, L2 



-ipDi/a \ r _ ( '^^Di - 1 -il){l - 2ai)vDi 
Di y ' ^ ~ V (1 - 2ai)vDi 2aDi 

2avDj — (1 — (7j)v —ipqDi — a 



qDj + il)j 2avDj + vaj 

where a = ip/ipi and q = {l — 2ai)v'^ /2 + 2 j ipj dl. The formula for the recursion 
operator Li is valid for arbitrary nonlinearity function '0(-^)i while operators L2 
and L3 arise for those functions ipil) that fulfill the condition Ref. [11] : 

^\ =0. (3.18) 

^iJ II 

It has been just this requirement which has defined our particular choice of the 
function il^{I) in initial statement of the problem. For this form of V' we have: 

<f = Ie-''\ a = b~\ q = v^/2-2ae-''^/b'^. 

In this particular case, the recursion operators read: 

^ _ f -e^^Di/a \ T - ( 2^//^ - 1 -e^^vDi/a 
^^ - \Di y ' ' " V ^Di 2Di/h 

( 2vDi/b-v 

^3 " V (^V2 - 2ae-''/b^) Dj + e-'^a/b 

^ebi (^2/2 _ 2ae-^^/b'^) Di/a - 1/6 
2vDi/b 

Let us now proceed with constructing the Lie-point symmetry group admitted 
by Eqs. (j3.1ip on the basis of these operators. 

An evident solution of the determining equation (|3.14p is 

/o=^, 9l = X- (3.19) 

The action of three recursion operators Lj, i = 1, 2, 3 on the vector with coordi- 
nates given by Eq. (j3.19p in accordance with Eq. (j3.16p generates the symmetry 
group given by 

/i = -e XI lot , /2 = ^- - T 



b a 

, 2vTi e*^ /w2 2aQ-^^\ x 

9i=Ti, gl = VTi + 2xi/b, 
1 fv"^ 2ae-^^\ ae-^^T 2vxi 



(3.20) 



Admissible for arbitrary nonlinearity ip{I), the symmetries /q, g'g describe the 
dilatation of r and x? whilst // and gl generate translations along u-axis. 

As it was formulated in Refs. [1HI20J. an invariant solution to the boundary 
value problem, in particular the one given by Eqs. (j3.1ip . must be found from 
the constructed Lie-Backlund symmetries under the invariance conditions 

/ = 0, g = 0, (3.21) 

supplemented by the original Eqs (|3.1ip . In Eqs. (j3.2ip the functions / and g are 
arbitrary linear combinations of coordinates ff and gf of the group generators 
Eqs. ()3.20p and must be chosen to satisfy the boundary conditions what in the 
actual case provides a localized intensity distribution. 

Unfortunately, the Lie point group generators p.20p with s = and s = 1 are 
not sufficient in order to determine a linear superposition able to satisfy a smooth 
localized (symmetric in x) intensity distribution at the boundary. Therefore, we 
shall continue using the discussed approach with operators Lj given by Eqs. (j3.16p 
and the vectors with coordinates of Eq. (j3.20p in order to find the Lie-Backlund 
symmetries of the higher order with s > 1. However, since the further calculations 
are quite cumbersome, it is convenient first to find the symmetry coordinates at 
the boundary where they have the simplest form. Afterwards, we completely 
reconstruct only those that will be included into the chosen linear superposition. 

Thus, at the boundary t = 0, v = the recursion operators read 

/ -e'^Di/a^ _(2Dj/h-l \ 

T -( 2/}?Di-l/h\ 

Action of these operators on Eqs. (|3.19p gives 9 symmetry operators Xi whose 
coordinates //, gf are listed in the Tabled) 

Based on the operators presented in the Table [H one can construct a lin- 
ear superposition providing a localized intensity distribution. For instance, the 
equation 

(2 - e'''-^)xii + 6(1 - ^''-^)XI = (3.24) 

has a particular solution x = \/2e~^^+^ — 1. Resolving / as a function of X; we 
get a convex symmetric on x intensity distribution 

/o(x) = ^l-ln(^^)) , (3.25) 

which is represented by a red curve in Fig. [2j 

From the Table [1] one can see that Eq. ()3.24p corresponds to the following 
superposition of symmetry operators of the first and the second order 

-al + y5? + ^gl = 0, (3.26) 



Table 1: Coordinates of the symmetry operators at the boundary. 




fi 


9i 


z = 





X 


i = 1 


-e'^xi/a 





i = 2 





IXlIb 


i = 3 


2xi/b'-x/b 





i = 4 


-2e''xii/ab 





i = b 





-ie^'^Xz/a - e^^x///a 


i = 6 





2X///&' - X//^ 


i = l 





4X7//^' 


i = 8 


-2e'"xii/ab - e^^x//« 





i = 9 


Axii/b^ - ^xi/b^ + x/b 





i = W 


Axn/b^ - 1X1 Ib^ 





i = 11 





2X///&' + XI lb 


i = l2 





-Aae-^^Xii/b'^ + Aae'^^i/b^ - axe""/^^ 



which, evidently, shall be supplemented by the equation: 

0. 



-/I + y/7 + y/2 



(3.27) 



It is easy to see that in order to find an invariant solution satisfying Eqs. (j3.1ip 
and the boundary conditions of Eq. (j3.25p . we have to reconstruct a complete 
form of the symmetry coordinates /I, g^ and /^ , g^. Acting by the operator Li 
on the couple fl, g\, we get 

-^^ Till a , 
-bQ^\xn^bxi)loL. 

The similar procedure applied to the operator L2 and coordinates /2, g\ yields: 
/2 = (4a _ Q^^v^}?)Tiila}? - 4t//6 - e^^vxija - Ae^^vxii/ab + r , 
g'^ = AvTii/b - VTi - e^^hv^xi/(^ + (4a - e^^ v^h'^)xii/oib'^ ■ 

These equations together with Eq. (|3.20|) represent the list of symmetry operators 
required for construction of analytical solutions. 



/I 



(3.28) 



(3.29) 



3.2 Invariant solutions 



Let us now find the desired analytical solutions. Equations (|3.26p - (j3.27p with 
expressions substituted from Eq. (j3.28p , ()3.29p represent a system of partial dif- 
ferential equations 

„,2l2 W\ / 2l2 6/\ 

"^" ^Xn + b(l-e''-'-^-4^)xi = 0, (3.30) 



2bvTii + I 2 - e' 
2vbe^^ 



bl-l 



a 



-XII 



2a 
a 



2a 



X/+ 2 



M-l 



2a 



Tii-bTi = 0. (3.31) 



The first integral to the equation (j3.30p can be easily found: 

2bvTi +1^2- e''-' - ^— J xi + bx = J{v). (3.32) 

J{v) in the above formula should be found from the comparison with Eq. (j3.3ip . 
Differentiating Eq. (I3.32p with respect to v, taking Eqs. (I3.1ip into account and 
comparing obtained expression with Eq. ()3.3ip . one can see that J{v) should be 
a constant. In view of a symmetric initial intensity distribution with respect to 
X — >■ —X reflections we are bound to choose J = 0. Then, substituting tj = —Xv 
into Eq. (|3.32p . we arrive at the following first order partial differential equation 

/ 2v2 bi\ 
- 2bvxv + 2 - e''-' - — )xi + bx = 0, (3.33) 



2a 
which can be integrated with a standard technique 



Integration of Eq. (I3.33|) gives two first integrals, 

Ji = ^, j^ = zl(2e^-'''-l + —vA. (3.34) 

V V \ 2a J 

Here and in what follows one has to keep in mind that a negative value of v 
corresponds to the focusing beam for the positive values of x. 

Now we are in a position to find a particular solution x{v, I) for the Eqs. ()3.1ip 
satisfying the boundary conditions Eq. (j3.25p . Let us first notice that from the 
system of equations (13. lip a linear second order partial differential equation 

aXvv + {e^^Xi)i = 0, (3.35) 

can be derived. Based on the result obtained in Eqs. (j3.34p . one can search for 
the solution to Eq. ()3.35p based on the following Ansatz 

X' = -vQ\J^), (3.36) 

Substituting Eq. (I3.36P into Eq. (|3.35p after some calculations we get: 

Qma* - (3/4 = 0, (3.37) 



where // = arcsinh( J^/y^ft^e/a). Equation (j3.37p has an evident general solution 

Q = Cie-''/^ + C2e''/\ 

where Ci and C2 are constants which should be found from the boundary condi- 
tions. 

Taking the Eq. ()3.36p into account we obtain the expression for x'- 

X = i-vy/^ [Cie-'^/^ ^ c'se'^/^j ^ (g^gg) 

10 



where e^'"^ is to be found from the equatfon 
sinh // = 



a —1 
Now we can express e^ from Eq. p.39p 



2e 



1-6/ 






e'' = K/2 ± yjKyl+l 



(3.39) 



(3.40) 



where 



V o^e V 






Summarizing, the following solutions to the system of equations (13. lip is ob- 
tained: 



X 



2ei-^^ - 1 + — — I + W ( 2ei-''^ - 1 + — — | + . (3.41) 



2a 



2a 



a 



In order to find the second function t{I,v), we shall integrate the original 
equation (j3.1ip keeping the result (I3.4ip in mind. From Eqs. (I3.1ip we have 



1 



„w. 



r = — / dve XI ■ 
a J 



For the sake of convenience, let us introduce a new variable ^: 



cosh(0 = 7 ^" . 



Then 



XI 



b ^^-^4+i-bl 



A^/Q sinh(0 

The expression for r becomes: 



,,irK^-i)'^%^,,/ZiJg!!_i^./fM_i 



?o 



2a 



(3.42) 



(3.43) 



& 



where a = \/2exp((l — hI)/2)). Taking the boundary conditions into account, 
the final solution reads: 



r = ,/-ln[e(''^-l)/^/^^^ + A/e''^- 



i (x^ + i) A 



(3.44) 



11 




Figure 2: Beam profile at different propagation distances, a = 3, b = 1. Red 
curve - z = 0, black curve - z = 0.8. 

After direct substitution of Eqs. ()3.4ip . (j3.44p into Eqs. (j3.1ip and a tedious 
calculation it is possible to verify that the obtained functions xi^iv) and t{I,v) 
are indeed exact analytical solutions for the formulated boundary value problem. 
In Fig.[2]we plot the intensity beam distribution at different propagation distances 
calculated on the basis of the found solutions Eqs. (|3.4ip . (|3.44p . 

Let us now examine the obtained result a little more closely. Firstly let us 
find the total radius of the beam as a function of propagation distance z. For 
this goal, it is necessary to determine where the intensity implicitly given by 
Eqs. (jOT]) . (fOIj) intersects the surface of x = 0. Putting / = in Eq. (jOij) . 
we find x\i=o = ±\/2e — 1. Substituting this number into Eq. (|3.4ip one can 
deduce that v = and, consequently, x = ±-v/2e — 1. Thus, the phase gradient 
at the beam edge is equal to zero, the total radius of the beam is a constant 
and does not depend on the propagation length. By the numerical integration 
of the solutions Eqs. (j3.4ip . (j3.44p one can verify their consistence with energy 
conservation, i.e. J I{x, z)dx from x = \/2e — 1 to x = —\/2e — 1 is a constant. 

The fact that the total radius of the beam for the case under consideration 
remains constant and does not depend on z is a new one and completely different 
from all exact analytical results obtained so far. It was demonstrated earlier 
[5tll2] that for the Kerr nonlinearity the total beam radius decreases upon beam 
propagation. In the present case, the beam shape and peak intensity are thus the 
only parameters depending on the propagation distance. 

Evolution of the beam peak intensity, which for the symmetry reasons has 
to be situated on the beam axis, can be easily found from Eq. (I3.44p . Putting 
X = 0, f = we have 



Wt = ^^" [{l/V2)e('^-'y' + V(l/2)e^^-i-l) • (3.45) 
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Figure 3: On-axial intensity distribution versus the propagation distance. 



On-axial intensity distribution versus the propagation distance is presented in 
Fig. [3] for a = 0.001, b = 0.2. We see that the intensity monotonically increases 
and tends to infinity for z approaching a critical value denoted as a self-focusing 
position Zsf. Its exact value can be found from direct analysis of the Eq. ()3.45p . 
Considering Eq. (j3.45p in the limit I — )■ oo one obtains 



Zsi 




(3.46) 



A detailed investigation of the Eqs. (j3.4ip . (j3.44p shows that the solutions exhibit 
no singularities before this point. 

4 Approximate solutions to the eikonal equation in 
(1+1) dimensions with arbitrary refractive index 

Because it is not possible to construct an exact analytical solution for every 
desired form of the refractive index and the boundary conditions, let us now 
present here a possible way to obtain approximate analytical solutions. We shall 
start from the Eqs. (j2.8p . Under certain conditions (see e.g. [E]), a can be 
considered as a small parameter. Due to its smallness, we will search for an 
approximate symmetries group operator: 

X =pdr + qd^, 

with coordinates in form of a power series in a: 

p = ^aV, q = "^a'q\ i = 0,...,oo. (4.47) 

13 



In case of Eqs. (|2.8p . the determining equations read: 



{Di + aDi)p-—-Djq = 0, {D^ + aDi) q + aDip = , (4.48) 

where w = v/a and 

oo oo 



s=0 



s=0 



Substituting Eqs. (j4.47p into Eqs. (j4.48p we arrive at a following system of recur- 
rent differential equations: 



Dlq^ + (1 - ,5,,o) [Diq'-^ + Dj{p^)) = , 
Dip' + (1 - 5,fl)Dip^-' - -Di{q^) = . 



(4.49) 
(4.50) 



The system (I4.49p - (l4.50p can be solved sequentially starting from a given g^. 
Thus, integration of Eqs. (|4.49p - (|4.50p gives 



■^ U=o 

oo 

q' = {l-ko) I dw\Y,T^+^d^^^q'-^ 
'' U=o 

oo 

- 9/ + EW"'X^+^^'% 



+- 



q')+P\ 



(4.51) 



s=0 



p'-')+Q' 



Here P* and Q* are arbitrary functions of invariants 



/, xh and Tl = T!-wY,C;{d^{Ilp)/dnxT'"'\ ^>1> 

p=0 

and Cp is the binomial coefficient (see Ref. |12j). 
Let us now put 

gO = o, / = l + 2/xx//v'- (4.52) 

This choice of p^ corresponds to a special form of initial light beam 



h 



{if/ 1) dl = -x^ 
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For nonlinearities (p{I) relevant to particular physical situations, such form of 
initial intensity distribution is very similar to the Gaussian profile. It should 
be stressed that, in general, there are no restrictions on the initial intensity 
distribution / = Io{x), the present choice is made only for the sake of further 
simplicity. 

Substituting this function into equation for g^ we get 

gi = wDj{2Ixxi/v) + G' = '^^\~ ^^ Di{21xxi/v) + G^ = -2(rx)/ , (4.53) 

ixi 

where the term with f was included into a new arbitrary function G^ , which was 
later put equal to zero. 

Using this result, we can calculate p in a similar way. We get 

p' = --{r'j) ■ (4.54) 

Finally, up to the first order of a, the coordinates of the approximate symmetries 
group generator read 

q = -2a{Tx)i, p=l + 2Ixxi/^-a^(^] ■ (4.55) 

Taking Eqs. (|2.8p into account, we rewrite the desired symmetry group operator 
as follows 

R= (l- ar^i^^^ J dr - 2xd^ + 2ar9/ . (4.56) 

Integration of the Lie equations corresponding to the point symmetry operator 
()4.56p gives us the three first integrals whose particular form depends on the 
choice of c/? = (/j(/). We write these integrals by introducing the new function 
4>{I), such that <pi = if/ 1: 

Ji = Xi J2 = ar^cpi -<j), ^ = 11] + —= I —j=-j=== . (4.57) 

V« J y<Piy/(P + J2 

Now, the solution is a function of these first integrals fulfilling the boundary 
conditions. This means that values of Jj are taken from the boundary conditions 
when r — )• and x^ x' ■ Hence we rewrite (j4.57p in the following form 

x' = x, 4>{Ioix')) = -az'^l'^(j)i + (j), 

'^^H'"'^'' A^ i AM (4.58) 

X y a \ =f + (x — fzjya 



0/\/0 + 'h J \/07\/<A + J2 

In the particular case of Kerr nonlinearity ip = 1, we arrive at the solution 
previously obtained in Ref. [T^ . 
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Let us now compare the approximate solutions with the exact solution con- 
structed in the previous section. In order to satisfy the boundary condition 
Eq. (j3.25p . p^ is taken in the form 

Then the generator (I4.56P reduces to: 

R = (26e + aT'^)dr - 2x9«, + 2ar5/ , 
and yields an approximate solution: 

({x - vzf + l) e^^ = al'^z^ + 2e , 




2a 1 / [a 

— -fa; — vz) arctg l^^ — 

e V V 2e 



(4.59) 



On-axial intensity distribution is given by expression e''-'^ — 2e = a/^z^, which 
exhibits no singularities. The on-axial intensity monotonically increases upon 
propagation. However, the function v{x,t) can exhibit singularities which, due 
to the symmetry of the problem, are expected to be on the beam axis. Let us 
investigate this behavior more closely. In the vicinity of the beam axis, Vx can 
be approximated as 

-- arctg ^/z^-J (^l--VTr^^*H V^Jj • 

Then, the position of a singularity at the beam axis (what corresponds to the 
rays intersection) can be found from the system of equations: 

2a 1 / , / « \ ^ .2 r2 _ o„ , „w 




l-ZsfJ — -aictgUzsfJ — j=0, azifl' = -2e + e''\ (4.60) 

The numerical solution of Eqs. (|4.60p yields Zgf ~ 1.03 6y^e/2a. As one sees, 
this result is very similar to Eq. (j3.46p obtained from the exact solution. The 
similar tendency has been observed in Ref. [21j for the Kerr refractive index: an 
approximate solution provides a longer self-focusing distance in comparison to 
exact one. 



5 Approximate solution to the Schrodinger equations 
with arbitrary refractive index in (1+2) dimensions 

In this section we shall turn to the construction of approximate solutions for 
the NLSE (j2.ip in (1+2) dimensions in media with arbitrary nonlinearity. Such 
a mathematical model describes the propagation of a continuum wave beam in 
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cylindrical geometry and has great number of application to particular physical 
situations (see e.g. [6ll8ll9] ). Let us begin with equations (|2.4p . ()2.5p supplemented 
by the boundary condition 



v{0,x) = 0, l{0,x) = N{x), 



(5.61) 



which corresponds to a collimated beam with arbitrary initial intensity distribu- 
tion. 

As usual, we start from construction of the Lie-Backlund symmetry operator 
of the form 

X = Kdv + \di . 

The determining equations read: 

DzH + vDxK + KVi — OLipDxX — aipjIiX 

- P [BiX + Bi^DxX + Bi^DlX + Bi^DlX) = 0, 



D^A + IDxK. + vDxX + Xvi + KIl^ \ = , 

X X 



(5.62) 



where 



B = Dr 



'Dx[xDx^I) 

x^fl 



' ^s ^^ "n « ' ^s ^ Tj T 
OX^ OX^ 



(5.63) 



Dx = dx + "^{vs+id^, + Is+idiJ, 

s=0 

and Dz we present as Dz = D^ + Dl , where 

oo ^ 

D'z = dz-J2(D'xivvi)dv.+ 

oo 

Dl = Y,[D'Ac^Vh+(3B)]dv.. 



D'x^\lv)+Dl{'-^^ 



di. 



s=0 



Since a and /3 in Eqs. ()2.7p can be considered as small parameters, we will 
search for k and A in form of a series expansion in powers of a and /3 



i,j=0 i,j=0 



ij) 



(5.64) 



and restrict ourselves only to the first order corrections 

K = K° + Kl + 0(a^/3^a/3), X = X^ + X^ + 0{a^, /3^,aP) , (5.65) 
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where k'> = k('>'^\ A^ = A^^'O), k' = aAc(i'O) + /3k(0'1), and A^ = aA(i'O) + /JA^O'^). 
Let us now write down the determining equations keeping only the hnear terms 
with respect to a and /3. We get 



(5.66) 



MoK° = , MiA° + M2k;° = , 

- 13 [S,A° + Bj^D^X' + Bi^^Dlx'^ + i?/...2?^°] = , 
MiX^ + DlX^ + M2K^ = 0, 

where 

Mo = D° + wZ)^ + VI , 

Ml = L»° + ?;Z)^. + vi+ v/x , (5.67) 

Ma = ID^ + Ii+ I/x . 

Let us now following to Ref. [LTj put 

^0 = Id.,{v^) , A° = -Z?.(/ra) . (5.68) 

Evidently this choice satisfies the zero-order Eqs. (j5.66|) and the invariance con- 
ditions: K*^ = 0, A*^ = at the boundary. Then k^ can be found from the first of 
the first-order equations in Eq. ()5.66p that is rewritten as 

Mo L' + a^h + f3D, ( ^^^l^^^ ) ) = 0- (5.69) 

The solution of this equation is expressed in terms of invariants of the operator 
Mo, 

where 

^ = di^, Six) ^ a^Nix)) + /? (^( V^^y)x )x (5_7i) 

Inserting this result into the second of the first-order equations in Eq. (I5.66P we 
get the equation for the function q^ 

MiX' + ^{IzD.,S{x)) = 0. 

X 

It is easy to show by direct substitution that the formula above can be rewritten 
as 

Mo{xX^) + aAl^D,S{x)) = 0. 

This equation can be integrated in a same as (j5.69p . Then one gets 



xX^ = —D^ixIzS 
18 



Xr 



Finally, up to the first order in the small parameters, the Lie-Backlund symmetry 
operators in the canonical form read: 

^ = ^^^1 + D, (s{x) -a^- p^I^^^^\ , (5.72) 



A = w ( Ji + - ) + Ivx - z 



I{l-zv,)S^^+{h + L]s^ 



. (5.73) 
We notice that based on ()2.4p . (j2.5p Eq. (j5.72p can be rewritten as follows: 

Together with Eq. (j5.72p . the equation above lead to two relations: 

v = zS^, (5.74) 

Vz = {I- Vxz)S^ , (5.75) 

that have to be fulfilled in order to preserve the invariance requirement k = 0, 
A = 0. 

Now, keeping the relation between the canonical form of the symmetries op- 
erator and the point symmetries group operator [1] in mind we can write down 
the group symmetry operator: 

i? = (1 + z^S^^) d, + S^d^ + {zS^ + vz^S^^)dx 

Operator Eq. (j5.76p is similar to the one obtained previously in Ref. p^ for a 
collimated beam with the exception that the S{x) now contains an arbitrary 
function $. The generator (|5.76p yields a system of characteristic equations: 

dz dv dx dln{Ix) . , 

A,A, A, A.A. 

This system of equations can be easily integrated after taking into account Eq. (|5.74p . 
The second and third equations of (j5.77p give 

S + ^ = S{i,), (5.78) 

where /i corresponds to the value of x s-t the boundary. The third and the fourth 
of Eqs. (|5.77p yield another invariant Ix/S^, what can also be rewritten as 

I = Nif,)^^, (5.79) 

X D^2 

where N{x) is an initial intensity profile. 
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x = x{l + 2z^S^2) , S{fi) = SU) + ^^, (5.82) 



From the first and the second of Eq. (j5.77p we have Vz = S^/{1 + z^S^y.). 
Taking Eq. I5.75P and DxS{x) = (1 — zVx)S^ into account and using 2xS^ = S^2, 
we finally arrive at a relation between x and x- 

x = x(l + 2^25^2). (5.80) 

Summarizing, the solutions are presented by the equations: 

v{x,z) = ^^^, I(x,z)=N{f,)^^, (5.81) 

where x 'ind fi are defined as functions of z and x via relations 

2 

These solutions describe the evolution of a collimated continuous wave laser beam 
with arbitrary initial intensity distribution in media with arbitrary nonlinear 
response. In case of Kerr refractive index, these solutions were investigated in 
details in Refs. J171I18J . obtained solutions to more complicated forms of the 
refractive index. 

We shall consider refractive index of the form 

n = no + ^2/ — (JkI ■ (5.83) 

Here, the second term in the right hand side represents the usual Kerr response, 
and the last term is responsible to the multiphoton ionization of the media in 
case of sufficiently strong electric field; the K then corresponds to the number of 
photons required for a simultaneous ionization, and ctk to the ionization rate. If 
at the entry plane of the nonlinear media the Gaussian intensity distribution is 
fulfilled, one obtains 

S = a exp(-x') - ajexpi-Kx^)/K + /3(x' - 2), (5.84) 

where 7 = KgkIq ~ jn^. 

Let us write the first of Eq. (|5.8ip as follows 

^ (5.85) 



1 + 2^25^2 



where 5*^2 = — ae"'^ + a^Q~^^ + /3. Here f is a single-valued function of x 
if the function x = x(^) can be determined from the first of equations (|5.82p 
uniquely. In order to find the region of multivaluedness let us investigate the 
function 1" = x(l + Iz^S^/z). We find where its first and second derivatives with 
respect to x vanish. 

From equation 1^ = we have 

1 + Iz^Sr^iri) + \z^T]S^rH = « ' (5-86) 
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where rj = x^- Equation Y-^-^ = gives 

x(35^„(r?)+2r/5^^^(77))=0. (5.87) 

Evidently, this equation is fulfihed if x = 0, or the expression in the brackets is 
equal to zero. In the first case, we observe the singularity at the beam axis, in 
the second case, the singularity takes place at the point rj = rjcr which should 
be found numerically from equation 35,,^ (r/) + 2?75^^,,(r/) = for each particular 
form of S{r]). 

In the first case, when x = 0; expression ()5.85p becomes singular on the beam 
axis at the point 



Zsf = l/V2(a(l-7)-/3). (5.88) 

The beam collapse at the beam axis occurs if a > aj + /3. This result is very 
similar to the case of Kerr nonlinearity considered in Refs. |17lll8j . the difference 
only coming from the presence of the 7 term under the square root. If a < a'y + /3, 
the beam collapse does not take place on the beam axis. 

In the second case, if x 7^ 0, position of the singularity can be found form the 
magnitude of rjcr which gives us x = i^Afcr and 5*^(7/) |^=^^^. Starting form Eq. 
(I5.86p . we can find the coordinates of the singularity position: 



-1 



^{jri + ■^^rjTi) 



±2^ < ^^^ 



V=Vcr 



(5.89) 

r)=r)cr 



Let us now consider two particular choices of the parameters in the refractive 
index. Let a = 0.01, /3 = 0.001, 7 = 0.1, K = 6. At the beam axis x = 0, the 
formula ()5.88p gives us the self-focusing position z^i ~ 7.9. Eq. ()5.87p has one 
solution 7] ~ 1.5, but the corresponding magnitude of the self- focusing distance 
defined by Eqs. (j5.89p is imaginary. This means that one observes only one self- 
focusing position at the beam axis. The beam intensity, I{x), and the phase 
gradient, v{x), as a function of x at different propagation distances are presented 
on Fig. m 

Let us now assume that a = 0.01, P = 0.001, 7 = 0.6, K = 8, then Eq. 
()5.87p gives us two solutions tji = 1.5 and r/2 — 0.11. Similar to the previous 
case, the first value of rj = iji gives no singularity. However, substituting r/2 into 
Eqs. (|5.89p . we get position where the beam collapse take place x ~ 0.1, z ~ 8. 
Considering behavior of the solutions at the beam axis, we see that the intensity 
increases, and the solution becomes singular at the point Zcr — 13, that already 
behind the distance at which the first singularity appeared. The intensity and 
phase gradient at different propagation distances are presented on the Fig. [5l 
Summarizing this part of the work, we can notice that for the form of the re- 
fractive index Eq. (|5.83p considered as an example in the present paper, several 
pictures in the global behavior of the solutions can be distinguished: i) if the so- 
lution of Eq. (I5.87P provide us with only imaginary values of Zgf j there is no beam 
collapse at all, ii) the singularity appears at the beam axis, ii) the singularity 
appears around the beam axis at the circle with radius given by Eqs. (|5.89p . 

21 





0.0 0.5 



1.0 1.5 2.0 2.5 



Figure 4: The beam intensity (left panel) and the phase gradient (right panel) 
distributions upon the transverse coordinate x for a = 0.01, (3 = 0.001, 7 = 0.1 
and K = 6. The transition near the beam axis a; — )• from the bottom to the top 
on the left panel and from the top to the bottom on the right panel correspond 
to curves at different distances z = 2, z = A, z = Q, z = 7.5. 
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Figure 5: The beam intensity (left panel) and the phase gradient (right panel) 
distributions upon the transverse coordinate x for a = 0.01, /3 = 0.001, 7 = 0.6 
and K = 8. The transition near the beam axis a; — )• from the bottom to the top 
on the left panel and from the top to the bottom on the right panel correspond 
to curves at different distances z = 2, z = 4, z = 6, z = 7.5. 

6 Conclusion 

In the presented paper, making use of the Lie symmetry analysis, we constructed 
exact and approximate analytical solutions for the problem of light propagation 
in highly nonlinear media. For the first time exact analytical solutions (j3.4ip . 
(I3.44P to the eikonal equations in (1+1) dimensions were found with nonlinear 
refractive index being a saturated function of intensity. It was shown that at a 
certain point at the beam axis Eq. (J3.46P the solution becomes singular: intensity 
tends to infinity asymptotically when the light propagation distance approaches 
Zsf = by^e/2a. 

In case of the eikonal equations with arbitrary nonlinear refractive index we 
constructed approximate analytical solutions. For the case of initial intensity 
distribution given by equation ()3.25p . the approximate solution was compared 
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with the exact one (j3.4ip , (j3.44p . It was shown that a value of the self- focusing 
position provided by an approximate solution was very close to the magnitude 
obtained from the exact formula (|3.46p . 

In the last section we considered a nonlinear Schrodinger equation in (1+2) 
dimensions with arbitrary refractive index. An approximate symmetry group 
admitted by both this equation and boundary conditions corresponding to col- 
limated beam with arbitrary initial intensity distribution was constructed. The 
solution was presented in the form of algebraic equations (j5.8ip which must be 
analyzed for each particular form of the refractive index and the initial inten- 
sity distribution. As an example, the case of two-term nonlinear refractive index 
Eq. (I5.83P was examined in details. We obtained an explicit formula for the self- 
focusing position Eq. (j5.88p and demonstrated that, for this form of the refractive 
index, the beam collapse can also take place outside the beam axis. 

Acknowledgment 

V.F.K and L.L.T acknowledge financial support by LiMat project. V.F.K also 
thanks DAAD for financial support. This work was also partially supported by 
RFBR projects No. 09-01-00610a and 08-01-00291a. 



23 



Bibliography 



[1] CRC Hanbook of Lie Group Analysis of Differential Equations, ed. 
N.H. Ibragimov, (CRC Press Inc., Boca Raton, 1994-1996). 

[2] V. F. Kovalev and D. V. Shirkov, Renorm-group symmetry for functionals 
of boundary value problem solutions, J. Phys. A: Math. Gen. 39, 8061-8073 
(2006). 

[3] V. F. Kovalev and D. V. Shirkov, Renormgroup symmetries for solutions 
of nonlinear boundary value problems, Physics Uspekhi, 51, 815 (2008). 

[4] Y.R. Shen, The principles of nonlinear optics, John Wiley &: Sons, Inc., 
New York - Chicester -Brisbane Toronto - Singapore, 1984. 

[5] S. A. Akhmanov, A. P. Sukhorikov, R. V. Khokhlov, Self-focusing and 
diffraction of light in a nonlinear medium, Sov. Phys. Usp. 10, 609-636 
(1968). 

[6] R. W. Boyd, Nonlinear Optics, (Academic Press, Amsterdam, Tokyo, 
2003). 

[7] J.-C. Diels, W. Rudolph, Ultrashort Laser Pulse Phenomena, (Elsevier, 
Amsterdam, 2006). 

[8] A. Couairon, A. Mysyrowicz, Femtosecond filamentation in transparent 
media, Phys. Rep. 441, 47-189 (2007). 

[9] L. Berge, S. Skupin, R. Nuter, J. Kasparian, J. -P. Wolf, Ultrashort filaments 
of light in weakly-ionized, optically-transparent media. Rep. Prog. Phys. 
70, 1633-1684 (2007). 

[10] S. A. Akhmanov, R. V. Khokhlov, A. P. Sukhorukov, On the self-focusing 
and self-chanelling of intense laser beams in nonlinear medium, Sov. Phys. 
JETP 23, 1025-1033 (1966). 

[11] V. F. Kovalev, and V. V. Pustovalov, Group and renormgroup symmetry 
of a simple model for nonlinear phenomena in optics, gas dynamics and 
plasma theory, Mathem. Comp. Modelling 25, 165-179 (1997). 

24 



[12] V. F. Kovalev, Renormgroup symmetries in problems of nonlinear geomet- 
rical optics, Theor. Math. Phys. Ill, 686-702 (1997). 

[13] L. Gagnon and P. Winternitz, Exact solutions of the cubic and quintic 
nonlinear Schrdinger equation for a cylindrical geometry, Phys. Rev. A 39, 
296306 (1998). 

[14] L. Gagnon and P. Winternitz, Lie symmetries of a generalised nonlinear 
Schrodinger equation: I. The symmetry group and its subgroups, J. Phys. 
A: Math. Gen. 21, 1493-1511 (1988). 

[15] L. Gagnon and P. Winternitz, Lie symmetries of a generalised non-linear 
Schrodinger equation. II. Exact solutions. J. Phys. A: Math. Gen. 22, 469- 
497 (1989). 

[16] V. I. Pulov, I. M. Uzunov and E. J. Chacarov, Solutions and laws of con- 
servation for coupled nonlinear Schrdinger equations: Lie group analysis, 
Phys. Rev. E 57, 34683477 (1998). 

[17] V. F. Kovalev, Renonrmalization group analysis for singularities in the wave 
beam self-focusing problem, Theor. Math. Phys. 113, 719-730 (1990). 

[18] V. F. Kovalev, V. Yu. Bychenkov, V. T. Tikhonchuk, Renormalization- 
group approach to the problem of light-beam self-focusing, Phys. Rev. A 
61, 033809(1-10) (2000). 

[19] L. L. Tatarinova, M. E. Garcia, Exact solutions of the eikonal equations 
describing self-focusing in highly nonlinear geometrical optics, Phys. Rev. 
A 78, 021806(R)(l-4) (2008). 

[20] D. V. Shirkov, V. F. Kovalev, The Bogoliubov renormalization group 
and solution symmetry in mathematical physics, Phys. Rep. 352, 219-249 
(2001). 

[21] V. F. Kovalev, Approximate transformation groups and renormgroup sym- 
metries, Nonlinear Dynamics, 22, 73-83 (2000). 



25 



